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ABSTRACT 

We present a general method to analyze reverberation (or echo) mapping data, that simultaneously 
provides estimates for the black hole mass and for the geometry and dynamics of the broad line region 
(BLR) in active galactic nuclei ( AGN) . While previous methods yield a typical scale size of the broad 
line region or a reconstruction of the transfer function, our method directly infers the spatial and 
velocity distribution of the BLR from the data, from which a transfer function can be easily derived. 
Previous echo mapping analysis requires an independent estimate of a scaling factor known as the 
virial coefficient to infer the mass of the black hole, but this is not needed in our more direct approach. 
We use the formalism of Bayesian probability theory and implement a Markov Chain Monte Carlo 
algorithm to obtain estimates and uncertainties for the parameters of our BLR models. Fitting of 
models to the data requires knowledge of the continuum flux at all times, not just the measured times. 
We use Gaussian Processes to interpolate and extrapolate the continuum light curve data in a fully 
consistent probabilistic manner, taking the associated errors into account. We illustrate our method 
using simple models of BLR geometry and dynamics and show that we can recover the parameter 
values of our test systems with realistic uncertainties that depend upon the variability of the AGN 
and the quality of the reverberation mapping observing campaign. With a geometry model we can 
recover the mean radius of the BLR to within ~ 0.1 dex random uncertainty for simulated data with 
an integrated line flux uncertainty of 1.5%, while with a dynamical model we can recover the black 
hole mass and the mean radius to within ~ 0.05 dex random uncertainty, for simulated data with a 
line profile average signal to noise ratio of 4 per spectral pixel. These uncertainties do not include 
modeling errors, which are likely to be present in the analysis of real data, and should therefore be 
considered as lower limits to the accuracy of the method. 

Subject headings: galaxiesractive — methods: data analysis — methods: statistical 



1. INTRODUCTION 

The energy emitted by active galactic nuclei (AGN) 
is argued to be the result of matter accreting onto 
supermassive black holes at the center of galaxies 
(|Lvnden-Bell fc Rees!ll971l ). However, the details of the 
geometry and kinematics of the region around the ac- 
cretion disk are not well understood . In the standard 
model of AGN, the region around the accretion disk is 
called the broad line region (BLR) due to the broad emis- 
sion lines from rapidly moving clouds of material nea r 
the black hole (|Antonuccilll993l ilUrrv fe Padovanilll995D . 
Models for the BLR attempt to explain the many cate- 
gories of AGN by the observer's viewing angle and the 
cov ering fraction of inflowing or o utflowing material (see 
e.g. iMurrav et aTlll995l; lElvisll2000t ). since the degree to 
which the BLR geometry and kinematics vary between 
individual systems is also unknown. 

The mass of the central black hole is a fundamental 
parameter in galaxy evolution, as suggested by a rela- 
tion between black hole mass and stellar velocity dis- 
persion of t he host galaxy bulge , the Mbh — rela- 
tion (see e.g. lBennert et al]|2011l ). While the black hole 
masses of very nearby galaxies can be measured by ob- 
serving the orbits of stars or gas, a different approach 
is needed for more distant galaxies because the gravi- 
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tational sphere of influ ence of the black hole c annot be 
spatially resolved (e.g. iFerrarese fc Fordl l2005t ). In ac- 
tive galaxies, the small size of the BLR, est i mated to be 
around ~ 10 14 - 10 16 m (jWandel et ai] |1999t iKaspi et all 
120001 : IBentz et al]l2006t ). inhibits direct imaging of the 
accretion disk and orbiting BLR clouds. 

Reverberation mapping provides a method to deter- 
mine the black hole mass, al ong with the geometry and 
the kinematics of the BLR (iBlandford fc McKedll982t 
iPetersonl 119931 : iPeterson et all 120041 ). Without reiving 
on spatially resolving the gravitational sphere of influ- 
ence of the black hole, reverberation mapping provides 
a powerful tool for studying black holes over a range of 
redshifts and black hole mass es (e.g . IPeterson et al.ll2004t 
IWoo et all 120071: IBentz et alJl2009t fDennev et al.ll2009f l. 
The method relies on the large time-variability of AGN 
luminosity, spanning t imescales of days to years (e.g. 
IWebb fc Malkanl [2000). Reverberation mapping data 
consists of a timeseries of frequent measurements of the 
intensity of the continuum and broad line emission. The 
line emission strength is assumed to be proportional to 
the continuum emission strength, but with a time lag due 
to the light travel time from the central ionizing source to 
the BLR material. The lag time between line and con- 
tinuum flux contains information about the size of the 
BLR, while the shape of the spectral line encodes the 
velocity information. An estimate of the black hole mass 
can be calculated assuming the BLR clouds orbit in the 
Keplerian potential of the black hole with velocities de- 
termined by the width of the spectral line and at a radius 
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given by the average lag between the line and continuum 
fluxes. 

A weakness of this traditional reverberation mapping 
method is that the relation between velocity and posi- 
tion observables of the clouds and the black hole mass 
depends on an unknown dimensionless proportionality 
constant that depends on the geometry and kinemat- 
ics of the BLR. In practice, this so-called "virial" fac- 
tor is estimated based upon some external criteria, such 
as the average factor that makes the Mbh — cr* rela- 
tion consistent between different black hole mass estima- 
tors and between samples of active and inactive galaxies 
(lOnken et al l 120041: iCollin et a.1.1 [2006L IW00 eFall l2l)10l: 
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2010l iGraham et al.l l201QfT Ideally, we 
would like to infer directly the morphology and kinemat- 
ics of the BLR, and the black hole mass, including its un- 
certainty (for a discussion of potential sys tematic errors 
in reverberation mapping see )Krolik|[200l . The models 
for the structure and kinematics may include a net inflow 
or outflow of BLR c louds, among other physically mo- 
tivate d models (e.g. iMurrav et al.l 119951 : iMarconi et al.l 
2008). Developing such a method is the goal of this pa- 
per. 

The data required for reverberation mapping encode 
the geometry and kinematic information to some degree, 
depending upon the quality, in the form of the transfer 
function (or response function) that maps the contin- 
uum emission onto the line emission. The average lag 
used to estimate black hole mass is the first moment of 
the transfer function. Previous analysis involved esti- 
mating the transfer function and then interpreting the 
transfer function in relation to a model of t he BLR (see 
Krolik fe Don! Il99l jPone fc Krolikl I199H: IHorne et al.1 



2003^ IBentz etall l201(l . This is necessary because the 



transfer function is a function of time lag, not position 
within the BLR, so its interpretation requires a physi- 
cal model. Estimating the transfer function requires in- 
ve rting a linear in t egral equa tion, and while t he me thod 
of iKrolik fe Donel (fl995h and lDone fc Krolikl (fl996h uses 
regularized linear inversion, thus allowing for uncertainty 
estimation, other inversion methods such as "maximum 
entropy" do not allow for straightf orward uncertainty 
estimation or mo del selection (e.g. IHorne et al.l 120031 : 
IBentz et al.ll2010l ). 

Our method of analyzing reverberation mapping data 
simplifies the process of obtaining a transfer function and 
then interpreting the result using different models. We 
compare reverberation mapping data directly with mod- 
els of the broad line region, obtaining uncertainty esti- 
mates as well as allowing for model selection. Once we 
have found models and model parameters that fit the 
data, we can easily compute the transfer function and 
average time-lag. Our goal is to constrain the geome- 
try and kinematics of the BLR and provide an internally 
consistent factor for the black hole mass. We note that 
the traditionally determined average time-lag is exactly 
equivalent to a model where the BLR is a face-on ring 
of a given radius (response = <5-function) or a spherical 
shell (response = step function). This implicit assump- 
tion drives the inference on the average lag and its result, 
as we will show in this paper. 

An important part of our method for directly model- 
ing re- verberation mapping data requires that we pre- 
dict the AGN continuum light curve between the ob- 



servations. Recent work has found that AGN contin- 
uum l ight curves are we ll-modeled by a damped random 
walk: iKellv et al.l (|2009[) used a con tinuous time stochas- 
tic process]_ Kozlowski et al.l (|2010[) used the formalism 
of lPress et all (|1992l ) and lRvbicki fc Press! (fl99l Il994| ). 
Th is model for AGN varia bility was applied to 900 AGNs 
by iMacLeod et al.l (|2010D in order to c orrelate v a riabil - 
ity with other parameters of AGNs. IZu et all ()2010f ) 
were then able to use this model for AGN variability 
to improve the standard analysis of reverberation map- 
ping data, including a better understanding of the uncer- 
tainties involved. They model the continuum light curve 
using Gaussian Processes to recover t he transfer func- 
tion, a ssumed to be a top-hat. As with iKozlowski et al.l 
(2010), they use an exponential covariance matrix to re- 
late the continuum flux at different points in the time 
series. We also use Gaussian Processes to model the 
continuum light curve, as well as a slightly more general 
exponential covarianc e matrix. Our method improves 
upon the approach of IZu et"ahl (|2010l ) by modeling the 
reverberation mapping data directly in terms of a geo- 
metric and dynamical model, rather than recovering the 
tr ansfer function . 

Bottorff et al.l ()1997l ) have also modeled reverberation 
mapping directly in an attempt to understand the BLR 
dynamics in the well-studied AGN NGC 5548. They ex- 
p and upon the hyd r omag netically driven outflow model 
of lEmmering et all ()1992f ) and use one set of parameter 
values to compare their model with NGC 5548. While 
the specific models presented here are clearly not as so- 
phisticated from a physical point of view, our method 
improves upon that approach by finding the best fit pa- 
rameter values of our simple models and believable esti- 
mates of their uncertainties. 

We consider two types of reverberation mapping data 
sets: velocity-unresolved, where there is a time series of 
the continuum flux and a time series of the integrated 
line flux, and velocity-resolved, where the data consist of 
a continuum flux time series, and a series of entire line 
spectra as a function of time. 

The paper is organized as follows. In §[2] we define 
and describe the physical problem. In §[3] we outline our 
methods in the formalism of Bayesian probability theory 
and describe the algorithms we use to compare reverber- 
ation mapping data to mock data created from a model 
of the BLR. In §|4]we test our method using simple mod- 
els of the BLR and show that we are able to recover the 
parameter values of our test systems. Finally, in §[5j we 
summarize our conclusions. Flux units throughout the 
paper are arbitrary, but computed consistently within 
our method. 

2. THE PHYSICAL PICTURE 

Throughout this paper, we assume a simple model for 
the BLR, described as follows. The AGN is defined to 
be at (0, 0, 0), and the observer is at (d, 0, 0). We model 
the distribution of BLR gas by defining the gas density 
profile p(x, y, z), assumed to be normalized such that 



p{x,y,z)dV = 1 



(1) 



where dV — dx dy dz and V is all of space. We assume 
that the gas absorbs the ionizing radiation, but is not 
self-shielding, so that gas at larger radii is still illumi- 
nated. It should be noted that our approach is fully 
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Fig. 1. — BLR clouds around the central ionizing source (central engine). The extra path length the light must travel from the central 
engine to the BLR cloud and then to the observer is the cause of the delayed response of the line flux. 




Fig. 2. — Simulated continuum emission datapoints with examples of the continuum interpolated using gaussian processes. The dispersion 
of the lines represents the uncertainty of the recovered light curve. As expected the unccrtanty is greatest where there are no data points. 
The top panel shows the simulated data used throughout this paper, whereas the bottom panel shows an example with gaps in the data. Our 
procedure takes into account the amount of information available and therefore the recovered light curve suffers from a larger uncertainty 
during the gaps. 



general and can support more complex models of the op- 
tical properties of the BLR, as well as its geometry and 
dynamics. 

2.1. Velocity- Unresolved Reverberation Mapping 

If the continuum flux varies with time according to 
/cont(i), then the total line flux as a function of time is 
given by 

/line W = A I /cont(i - l(x, y, z))p{x, y, z) dV (2) 
Jv 



where l{x,y,z) is the lag, or time delay, associated with 
BLR gas at position (x,y,z), and A is a response coef- 
ficient. The lag I for each position is simply the excess 
light travel time from taking a path starting at (0,0,0) 
that travels to some gas at (x,y,z), where the light is 
absorbed and reemitted as line emission, and that finally 
travels to the observer, relative to a direct path straight 
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from the AGN to the observer: 

l(x, y, z) = (v 'x 2 + y 2 + z 2 (3) 

+ y/(x - d) 2 + y 2 + z 2 - dj /c (4) 



For any case of interest, d x 2 + y 2 + z 2 , and there- 
fore this is well approximated by: 



l(x, y, z) w (\/ x 2 + y 2 + z 2 — x\ jc 



(5) 



which is the formula adopted throughout this paper. See 
Figure Q] for an illustration of this model. 

Note that Equation [5] is a special case of the general 
equation 



/line(t) = A J *(r)/ cont (i - r) dr 



(6) 



where ^(t) is the so-called transfer function, which gives 
the response of the line flux to a delta-function pulse in 
the continuum flux 3 . Thus, for any particular system, if 
we can infer the density of BLR clouds throughout space, 
we can automatically deduce the corresponding transfer 
function: 



5(t- l(x, y, z)) p(x, y, z) dV 



(7) 



The meaning of this equation is that each location in 
space contributes to the transfer function at the value of 
the location's lag, with the size of the contribution being 
proportional to the amount of gas at that location. 

2.2. Velocity- Resolved Reverberation Mapping 

Now suppose that the BLR gas is in motion, such that 
the system can be described by a time-invariant distri- 
bution function g defined over the phase space of a single 
particle: 

g(x,y,z,v x ,v y ,v z ) = p(x,y,z)g(v x ,v y ,v z \x,y,z) (8) 

The motion of the gas along the line of sight is assumed 
to affect the wavelength of reemitted light, but its dis- 
tribution function is assumed to be time invariant and 
therefore does not vary during the observing campaign. 
Then the emission line profile at time t will be a function 
of the line of sight velocity, v\ os : 



f\in e (vios,t) = A / / / co nt(* ~ l(x,y,z)) (9) 

J Vy : V Z J V 

xg(x, y, z, v x , v y ,v z ) dx dy dz dv y dv z (10) 

where v\ os is in the x direction. This is the velocity- 
resolved equivalent of Equation [2] 

3. METHOD 

Our method for constraining the geometry and kine- 
matics of the BLR is an application of Bayesian Inference 
(jSivia &: Skillingl l2006f ) . In general, to infer parameters 
8 from data D, we begin by assigning a prior proba- 
bility distribution p(8) describing our initial uncertainty 
about the parameters. Sampling distributions p(D\8) are 
also assigned to describe our uncertainty about how the 

3 For readers more familiar with image analysis, the transfer 
function is analogous to a PSF. 



data are related to the parameters. Once specific data 
D = D* are obtained, our updated state of knowledge 
about the parameters is described by the posterior dis- 
tribution, given by Bayes' rule: 



p(6\D = D*,I)<xp(8\I)p(D\6,I)\ 



D=D' 



(11) 



Here I is any background information we have about 
the problem. In complex problems, where 8 consists 
of a large number of parameters, Monte Carlo methods 
are used to produce random samples from the posterior 
distribution for 8. Methods suc h as Nested Sampling 
(|Brewer. Partav, &; Csanvil l2010f) can also provide the 
normalization constant for the posterior, known as the 
evidence, which is the key quan tity for comparing th e 
entire model with an alternative ([Sivia fc Skilling 2006). 

In our method, the parameters 8 to be inferred are 
those describing the spatial profile of the BLR gas, and 
the continuous continuum flux, / CO nt(£)- Since it is im- 
possible to represent a continuum in a computer, we in- 
stead infer / con t(i) evaluated at 500 time points, covering 
a time interval larger than the continuum data. The con- 
tinuum modeling technique is described in detail in the 
next section. 

Throughout this paper, both the continuum flux and 
line flux timeseries are considered part of the dataset D: 



D {yiine 5 ycontinuum} 



(12) 



The prior information consists of the times at which the 
line flux and continuum flux are measured, t and the 
error bars on the line flux and continuum flux measure- 
ments, <t: 



^ = {(t J 0-)l in e.( t '°-) 



/continuum 



} 



(13) 



The likelihood function is chosen to be Gaussian, cen- 
tered around the model-predicted line flux timeseries: 



p(D\8) 



n exp 
II 



yi,Unc-mi(8) 



(Ki7i)v27r 



(14) 



where n is a "noise boost" parameter to account for the 
presence of unknown systematic effects not included in 
the reported error bars, such as those due to flux calibra- 
tion, wavelength calibration, and continuum subtraction. 

Once the posterior distribution is obtained, many dif- 
ferent algorithms are available for exploring it and com- 
puting summaries such as marginal distributions for pa- 
rameters. We have implemented our model with two 
methods, the first is Metropolis-Hastings, a Markov- 
Chain Monte Carlo (MCMC) algorithm, which pro- 
vides samples from the posterior PDF for the model 
arameters. The second is Diff usive Nested Sampling 
Brewer. Partav. fc Csanvil [20101. which provides sam- 
ples from the posterior PDF and an estimate of the 
evidence value for the model. Although the evidence 
calculation makes the second algorithm significantly 
slower than the first, Diffusive Nested Sampling is much 
faster than altern ative MCMC-based implementat ions of 
Nested Sampling (| Brewer. Partav. fc Csanvil [20101. The 
results presented here to test the method use the MCMC 
algorithm, while the Diffusive Nested Sampling algo- 
rithm is used t o apply the method to real reve rberation 
mapping data ([Brewer. Treu. fc Pancoastl|2011l in prep). 
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3.1. Continuum Interpolation 

In order to create a mock line flux time series to com- 
pare with the data, it is necessary to interpolate between 
the continuum flux datapoints. Linear interpolation is 
the simplest approach, but it does not provide an esti- 
mate of the uncertainty in the interpolation, suggesting 
that we know precisely the value of the continuum f(t) 
at all times between the measured datapoints. If we want 
to obtain reliable uncertainties in our results, we should 
acknowledge the uncertainty introduced by the interpo- 
lation process. 

To account for this, we consider the entire contin- 
uum function f con t(t) to be an unknown parameter 
to be inferred from the data. The prior distribu- 
tion for fcont(t) is a Gaussi an Process ( M acKavl 120031 : 
iRasmussen fc Williams![2006T ) . which is a convenient class 
of probability distributions over function space. Given a 
mean function fi(t) and a covariance function C(ii,t2), 
the probability distribution for the function value / at 
any finite set of times is a multivariate Gaussian: 



p(f\n > C) = 



1 



v/(27r)"detC 



exp 



> T c-\f-n) 

( 15 ) 

where \x is a vector of means at the relevant time-points, 
and C is the covariance matrix, obtained by evaluating 
the covariance function at the relevant times. In the re- 
verberation mapping problem, f con t(t) is constrained by 
two data sets: the continuum measurements, and the line 
measurements. We parameterize the covariance function 
and mean function with four hyperparameters: /i (the 
long-term mean), a (the long-term standard deviation), 
r (typical timescale of variations) and a (a smoothness 
parameter between 1 and 2), such that the mean func- 
tion is a constant fj,(t) = fi and the covariance function 
is 

' 1*2 — *lT 



C(h,t 2 ) 



a 2 exp 



(16) 



The posterior distribution function for f(t) given some 
continuum data (but not the line data) is shown in Fig- 
ure [2j Note that outside the areas where we have data, 
the uncertainty gets large, but in areas where the data 
are well sampled, the uncertainty in the interpolation is 
small. We keep track of f(t) at 500 times, both slightly 
preceding and following the data. Further interpolation 
between these 500 points is linear. 500 continuum param- 
eters is sufficient to render the distance between contin- 
uum flux points much smaller than the maximum moni- 
toring cadence, allowing us to resort to linear interpola- 
tion only on scales not probed by the data. We change 
the 500 parameters in the same way as the model param- 
eters, with every new proposal for the continuum func- 
tion related to the one before. The function f(t) can be 
parameterized by 500 variables with standard normal pri- 
ors, which are converted to f(t) values by multiplication 
with the Cholesky decomposition of C. We note that our 
Gaussian Process method for interpolation, in the spe- 
cial ca se a = 1, is equivalent to the method of IZu et all 
( 2010), apart from computational details, a = 1 has also 
been used in detailed studies of quasar variability (e.g. 
IMacLeod et~aT1l2010D . 

3.2. Creating Mock and Simulated Data 



Given the phase-space density for the BLR gas and 
the continuous continuum light curve, we can easily cre- 
ate a mock line flux timeseries by adding together the 
line flux from all the gas, which is proportional to the 
continuum flux at the respective lag of the gas. The re- 
sulting mock line flux timeseries can then be compared 
to the reverberation mapping data and does not depend 
on the kinematics of the gas. If we include the velocity 
information of the gas, we can create a mock spectrum 
for each point in the timeseries. In order to create a 
mock spectrum, we make a histogram weighted on flux 
of the amount of gas with a given velocity using the same 
velocity resolution as the data. We then convolve the 
histogram with a gaussian whose width is defined by a 
combination of thermal broadening and instrumental res- 
olution. The mock spectrum can then be compared to 
the reverberation mapping spectral data and depends on 
the kinematics of the gas. 

4. ILLUSTRATION AND TESTS USING SIMPLE MODELS 

In order to illustrate our method, we have developed 
simple models of the BLR region geometry and dynam- 
ics. As this method is fully general, it is also possible 
to implement more complex models within the frame- 
work described so far. We showcase these simple models 
by creating simulated data with known true parameter 
values in our models. This allows us to test our code 
as well as to explore the accuracy and precision of the 
results obtainable by this method of reverberation map- 
ping analysis. Such tests on simulated data also allow 
us to ascertain the data quality needed to perform infer- 
ences regarding increasingly complicated model param- 
eters. We showcase both geometry-only and geometry 
plus kinematics models, where the latter are the same as 
the first with the addition of velocity information given 
to the BLR gas. We show transfer functions for the ge- 
ometry models and velocity-resolved transfer functions 
for the kinematics models. 

To ensure that in our method the true parameter values 
are recovered, we save instances of each model and use 
them as simulated reverberation mapping data, adding 
noise and varying the timeseries characteristics to match 
reverberation mapping campaigns of varying quality. A 
simulated dataset consists of line flux and continuum flux 
measurements. Given a BLR model, the continuous con- 
tinuum light curve is all that we need to create, since the 
mock line flux measurements can be obtained from the 
model and continuous continuum light curve. We create 
continuous continuum light curves by using the hyperpa- 
rameters of the Gaussian Processes continuum interpo- 
lation. The hyperparameters contain information about 
the timescales and levels of variability in an AGN contin- 
uum timeseries. We use values for the hyperparameters 
from int erpolation of the Lick AGN Monitor ing Project 
(LAMP: rWalsh et al.ll200l iBentz et al.ll2009[ ) continuum 
timeseries of Arp 151, one of the most variable AGN in 
the LAMP sample. The values used for the hyperparam- 
eters were \i = 75 (arbitrary units), a = 30 (same units 
as /i), r = 6 x 10 6 seconds and a = 1.5 (dimensionless) . 

4.1. Geometry Model: Ring/Disk/Shell 

4.1.1. Model definition 

We use a flexible geometry model of the BLR gas den- 
sity to test our method when only integrated line flux 
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Fig. 3. — Example spatial distributions of the broad line emitting 
gas that can be recovered by our generic geometric model. They 
include a ring/disk (top panel), a spherical shell (middle panel), 
and a spherical gaussian distribution (bottom panel). 



measurements are used instead of the full spectral shape. 
The model is that of a spherical shell centered on the 
central engine with parameters allowing partial, axisym- 
metric illumination of the shell and varying inclination of 
the resulting ring/disk. Examples of possible configura- 
tions, ranging from a complete shell to a thin ring/disk, 
are shown in Figure |3] The parameters of the model are 
the mean radius of the disk, ro , the thickness of the disk 
in the radial direction, ay, the illumination angle of the 
shell, and the inclination of the shell. The illumination 
angle is defined so that values approaching define an 
increasingly thin ring/disk and a value of 7r/2 defines a 
spherical shell. The inclination angle is defined so that 
values approaching define a face-on ring/disk and a 
value of it/2 is an edge-on ring/disk. We use a normal 



distribution to define the radial thickness of the shell, so 
that ro and ay are the average and la width of a normal 
distribution. The normal distribution is created in the 
x, y, and z cartesian coordinates. 

It is important to set appropriate prior probability dis- 
tributions for each model parameter. For parameters 
where we know the order of magnitude of the parameter 
value we use a flat prior in the parameter. Examples of 
parameters with flat priors in the parameter include the 
inclination angle and the illumination angle, which may 
only vary between and it/2. For parameters where we 
do not know the order of magnitude of the parameter 
value we need a prior that treats many orders of magni- 
tude equally, so we use a flat prior in the log of the pa- 
rameter. Examples of parameters with flat priors in the 
log of the parameter include ro and ay. These choices of 
prior probability express complete ignorance in the value 
of a parameter within some reasonable range, but it is 
necessary to make the distinction between whether or not 
the order of magnitude of a parameter value is known. 
In the cases considered in the remainder of this paper, 
the posterior is much narrower than the prior, and there- 
fore the inference is dominated by the likelihood, i.e. the 
information contained in the data. 

The underlying spherical symmetry of these models 
and the angular dependence of the ring/disk model allow 
us to use spherical coordinates. In order to sample the 
gas phase-space density at a finite number of points, we 
use a grid in log(r), <j>, and cos(0). Using equal steps in 
cos(0) instead of 9 means that the volume of each grid 
point depends only on the radius, r. The density is then 
multiplied by the volume of the grid point to find the 
total mass of gas in each grid point. The emissivity of 
each grid point also depends on the radius r because the 
continuum ionizing radiation flux falls off as r~ 2 , requir- 
ing more gas mass at larger radii to have the same line 
flux contribution as less gas mass at smaller radii. In 
general, the illumination parameter allows us to model 
any axisymmetric ionizing flux. 

We test our method to recover the BLR model pa- 
rameters by creating simulated data, the true parame- 
ter values of which are given in Table [TJ The continous 
continuum function is obtained using the hyperparame- 
ters from the Gaussian Processes interpolation of Arp 
151 reverberation mapping data, as described in Sec- 
tion ^. 1[ and evaluated at 120 consecutive "observations" 
one day apart. The line flux timeseries for each model 
are generated using this continuous continuum function 
and a given set of model parameters. The line flux time- 
series contain 60 "observations" one day apart, starting 
60 days after the start of the continuum flux "obser- 
vations" . These simulated data are meant to represent 
excellent reverberation mapping data, with an observa- 
tion camgmgjiofsjmilar length to recent campaigns (see 
e.g. iBentz et al.ll2009t ). but without gaps due to difficult 
weather conditions. Additional noise has also been added 
to the simulated data. Most simulated datasets have 
line flux errors of 1.5%, which represents very favorable 
observing conditions, but we have also tested simulated 
data with errors of 5% to reflect the current typical error 
of reverberation mapping line flux measurements. 

4.1.2. Testing the geometry model 
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Fig. 4. — Posterior probability distributions for face-on disk ge- 
ometry model parameters of simulated data 4 (see Table [TJ with 
1.5% line flux uncertainty. Top to bottom: ro, ay, inclination an- 
gle, and illumination angle. The inclination angle and illumination 
angle both have a resolution given by the grid in cos 0. The true 
value for each parameter in this model is shown by the vertical 
red line and the grid in cos 8 is shown along the x-axis with green 
crosses for the angular parameters. The grid used to create these 
posterior distributions is 60 steps in log(r), 40 steps in <f>, and 60 
steps in cos 9. 



Fig. 5. — Posterior probability distributions for shell geometry 
model parameters of simulated data 5 (see Table [TJ with 1.5% 
line flux uncertainty. Top to bottom: ro, oy, inclination angle, and 
illumination angle. The true values for the parameters and the grid 
points are shown as in Figure[4] Note that since the simulated data 
is spherically symmetric, it should not strongly prefer an inclination 
angle, and thus no true parameter value is shown in the inclination 
angle pdf. 



The first test is whether we can recover the parameter 
values of the simulated data using the MCMC algorithm 
described in Section [3] Since our one flexible geome- 
try model encompasses a number of different geometries, 
such as a shell, thin or thick ring or disk, we do not 



have to consider model selection at this point. We test 
the many possible geometries of this model by creating 
five simulated datasets, whose true parameter values are 
given in Table [T] The simulated datasets include an in- 
clined disk with line flux errors of f .5% and 5% and an 
edge-on disk, a face-on disk, and a shell with line flux 
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Fig. 6. — Joint posterior probability distributions for inclination 
and illumination angles for face-on disk with 1.5% line flux uncer- 
tainty (simulated data 4) and shell with 1.5% line flux uncertainty 
(simulated data 5). The true parameter values are shown by (top) 
the black cross and (bottom) the black dashed line. 



Fig. 7. — Timeseries for face-on disk (simulated data 4, top 
panel) and shell (simulated data 5, bottom panel), both with 1.5% 
line flux uncertainty. Simulated data are shown in blue with error 
bars and the mock data from a random set of parameter values 
sampled from the posterior is shown in red. The continuum light 
curve used to create these line light curves is shown in Figure[2] 



errors of 1.5%. The MCMC algorithm is typically run 
for 150,000 iterations and all parameter values are recov- 
ered to within two standard deviations of the posterior 
probability distributions of the parameters, with 10/13 
recovered to within one standard deviation. This is as 
expected, since we should find the true parameter value 
to lie within la about ~ 68% of the time and to lie within 
2cr about ~ 95% of the time. The mean and standard de- 
viation of the posterior distributions are given in Tabled 
with the exception of many of the angular parameters, 
where the quoted la uncertainty does not adequately 
describe the posterior distribution. Part of the reason 
for the standard deviation of the angular parameters not 
describing the posterior is due to the uneven step size 
in 9, so that values of the illumination angle close to 
7r/2 and values of the inclination angle close to radians 
have much poorer angular resolution. This might lead 
to an angular parameter being quoted as having a mean 
of 1.22 radians and a la uncertainty of 0.29, as for the 
illumination angle of the Shell model simulated data, but 
while this uncertainty may seem large, it corresponds to 
an uncertainty of only 1-2 grid points in 9. The poste- 



rior distributions for the face-on disk and shell simulated 
data are shown in Figures [U and \5\ Select joint probabil- 
ity distributions between the inclination and illumination 
angles are also shown in Figure [6] in order to show the 
degeneracies between different models. In particular, for 
the shell model, the inclination is not constrained un- 
less the illumination angle is small, or rather, unless the 
sphere of BLR gas is not entirely illuminated. 

The posterior pdfs for the five simulated datasets show 
that the edge-on disk, face-on disk, and shell geometries 
allow for excellent recovery of the parameter values with 
estimates of the uncertainty. For the two inclined disk 
simulated datasets, there is some degeneracy in the an- 
gular parameters, leading to large uncertainties in their 
average values. The MCMC algorithm finds a more likely 
geometry configuration than the true configuration for 
the inclined disk datasets, although the true configura- 
tion is still a valid possibility with posterior local max- 
ima at the true parameter values. With the increased 
simulated line flux error from 1.5% to 5% however, it 
becomes increasingly difficult to recover the angular pa- 
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Fig. 8. — Velocity-unresolved transfer functions for face-on disk 
(simulated data 4) and shell (simulated data 5), both with 1.5% line 
flux uncertainty. The same grid was used to make these transfer 
functions as was used to obtain the posterior probability distribu- 
tions shown in Figures[4] and \E\ 



ramcters, and only the mean radius is recovered with a 
small enough uncertainty as to be useful in describing 
the BLR. This emphasizes the importance of obtaining 
high quality line flux data in reverberation mapping cam- 
paigns. 

The timeseries and transfer functions for the face-on 
disk and shell MCMC geometry model tests are shown 
in Figures [7] and [51 respectively. The timeseries show 
the simulated data overlaid with mock data created with 
parameters sampled randomly from the posterior proba- 
bility distributions. The fit of the mock data to the sim- 
ulated data is excellent for all five models. The variety 
in the shape of the simulated data timeseries, all well-fit 
by their respective models, shows that the MCMC algo- 
rithm for model parameter value recovery is robust for a 
wide range of models. The transfer functions also show 
a variety of shapes. For a thin shell geometry, thinner 
than the shell of simulated dataset 5, our resulting trans- 
fer function agrees with th e analytic form of a tophat 
function fsee lPetersonlll993l) . 

4.2. Dynamical Model 




Fig. 9. — Sketch of the dynamical model. The angular momen- 
tum vector L defines the plane of the orbits. Owing to cylindrical 
symmetry, for each value of 8g we consider the entire family of L 
generated by rotation around the z-axis. The observer is assumed 
to be in the x-z plane, at angle Oi from the z-axis. 
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Fig. 10. — Illustration of the combined constraints given by the 
illumination function and by the dynamical model. The red line 
shows an example of the distribution of illuminated BLR gas mass 
assuming a uniform underlying density. The blue line shows the 
actual underlying mass distribution as constrained by the dynami- 
cal model. The resulting effective distribution of illuminated mass, 
consistent with both the geometry and dynamical constraints is 
given by the product of the two functions, shown in black. 



4.2.1. Model definition 

In order to constrain the kinematics of the BLR and 
the mass of the central black hole, we must model the ve- 
locity distribution of the BLR gas in the context of a dy- 
namical model. For simplicity of illustration and speed of 
computation, we consider here a cylindrically symmetric 
model where the BLR gas is considered to be made of test 
particles in bound orbits within the spherical Keplerian 
potential of the black hole. We parameterize the model 
in terms of energy and angular momentum, constants of 
the BLR gas motion, so we are guaranteed velocity and 
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Fig. 11. — Example spectra from three simulated datasets: 
(top) face-on disk with orbits confined to the disk, (middle) face- 
on disk with isotropic distribution of orbit orientations, and (bot- 
tom) spherical distribution with isotropic distribution of orbit ori- 
entations. The instrumental resolution of the simulated spectra 
is FWHM ~ 800 kms -1 . The bottom spectrum for a spherical 
distribution of orbits is wider than for a face-on disk because the 
spherical distribution allows for orbits to move directly along the 
line of sight, while the face-on disk only results in a small com- 
ponent of the BLR gas velocity lying parallel to the line of sight. 
The width of the spectral line is thus directly connected to both 
the opening angle of the disk and the inclination angle. 



geometry distributions that do not evolve in time, and 
are therefore stationary during the monitoring campaign. 
In future papers we will generalize the model to include 
unbound orbits to describe inflows and outflows, and also 
other physical mechanisms, such as radiation pressur e or 
winds (jMarconi et al.ll2008i: iNetzer fc Marzianill2010[) . 

The model is illustrated in Figure [§] For any choice 
of angular momentum L, energy E and black hole mass 
Mbh, the motion of the BLR test particles is then de- 
scribed by the standard conservation equation resulting 
in elliptical orbits in the plane perpendicular to the an- 
gular momentum. Given our cylindrical symmetry we 
will consider families of angular momenta obtained by 



0.25 




Log[Mean width/m] 

Fig. 12. — Posterior pdfs for the first dynamical simulated 
dataset: face-on disk with the orbits confined to the disk. (Top) 
black hole mass, (middle) the average radius of the BLR gas mass, 
and (bottom) the average width of the BLR gas mass. 



rotation around the z-axis and defined by the polar an- 
gle #o (see Figure [9j . The spatial density of the BLR is 
then given by 

P(r, 6, <j>\E, L, »„)a-x 1 — (17) 

v |\/sin 2 6>o - cos 2 (9 1 

where the angular term comes from integrating over the 
uniform distribution of azimuthal angle 0o of the angular 
momentum vector, and v is the total magnitude of the 
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velocity vector: 



2E 



1GM BH 



(18) 



Owing to the symmetry of our model we can consider 
only 6q < n/2 (i.e. L z > 0), obtaining the following 
limits on the allowed 9 coordinate for the BLR: 



- do < 9 < 



(19) 



As 9q approaches zero, the model represents a thin disk, 
while as 6q approaches tt/2 the model covers the whole 
sphere. Conservation of energy and angular momentum 
limits the radial coordinate to the range: 



r > 



r < 




(20) 



(21) 



Finally, E and L are connected by the usual condition: 



\L\ < 



GM BH 



V^2E 



(22) 



For every allowed value of r, 8, and <f>, the component 
of the velocity vector along the line of sight can be com- 
puted in the standard manner, resulting in two solutions 
per position, in general (outbound and inbound; four if 
one considers L z < as well). More complex geome- 
tries and kinematics can be obtained by superpositions 
of multiple sets of E, L, and 9q values within the same 
potential given by Mbh- However, this further increases 
the dimensionality of parameter space and computing 
time. Therefore in this example we will only use one 
such set. 

We apply prior probability distributions to the model 
parameters as described for the geometry model. The 
priors for the extra parameters in the dynamical model 
not part of the geometry model are as follows. The pa- 
rameter #o has a flat prior in the parameter ranging from 
to 7r/2. The parameters Mbh, E, and L have flat priors 
in the log of the parameter. 

In addition, in order to impose a BLR gas geome- 
try, we model the distribution of illuminated gas, as the 
product of the spatial distribution given by the dynam- 
ical model with that imposed by one of our geometrical 
models, representing in this case the illumination func- 
tion. This results in a broad range of geometries, giv- 
ing the model a considerable flexibility (for example, in 
the future one could think of an anisotropic illumination 
function to model dust obscuration). The procedure is 
illustrated in Figure [10] Note that the radial distribution 
of the illuminated gas is not gaussian anymore, as in the 
ring/disk/shell geometry model. The mean radius then 
is not the ro parameter of the geometry model, but must 
be computed numerically for each set of geometric and 
dynamical parameters. Similarly, the mean width is no 
longer ay and must be computed numerically. 

A model spectrum at a given time is obtained by sum- 
ming all the line of sight velocities, weighted by the spa- 
tial density of illuminated gas multiplied by the contin- 
uum flux at an epoch corresponding to the appropriate 



lag-time. In order to compare with real data, the model 
spectrum is then convolved with a gaussian to represent 
instrumental broadening. Since we do not expect real 
data to match our model perfectly, we introduce a rela- 
tively large uncertainty in the form of the spectral line by 
adding gaussian noise with a variance of <J 2 {F) = a F+/3, 
where a = 0.00018 and /3 = 0.025. This model for the 
variance assumes both a dependence on spectral line flux 
F through the a parameter and a dependence on the con- 
tinuum uncertainty through the /3 parameter. The units 
of a are flux and the units of f3 are flux 2 . The specific 
values of a and are related to the arbitrary flux units 
of our simulated spectra and result in a signal to noise 
of ~ 4. Conservatively this signal-to-noise ratio is lower 
than typically achieved in state of the art spectral mon- 
itoring campaigns. Examples of synthetic spectra at a 
resolution of FWHM= 13. lA, or'- SOOkms" 1 at the 
wavelength of H/3, are shown in Figure [TTJ The face-on 
disk systems (top and middle panel of Figure [TT] have ve- 
locity bins of — 120 km s -1 while the sphere system (bot- 
tom panel) has velocity bins of — 20kms~ 1 . Notice how 
the line shapes are clearly different even for models with 
the same black hole mass. This is a clear illustration of 
the power of velocity resolved reverberation mapping as 
a diagnostic of the BLR geometry as well as kinematics. 

4.2.2. Testing the dynamical model 

We test our dynamical model by creating simulated 
data-sets consisting of timeseries of the continuum flux 
and of the line profiles of a broad line. The line profiles of 
the simulated datasets are shown in Figure[TlJ The kine- 
matics parameters E and L are initially chosen to satisfy 
nearly circular orbits of the BLR gas at the mean ra- 
dius given by the illumination function. A disk of broad 
line emitting material can be constrained by either the 
illumination function or the value of 9q . 

The first simulated dataset is a thin disk viewed nearly 
face-on, with dynamics imposed by a single value of 
energy and angular momentum. The thin disk is con- 
strained by the value of 9q, while the illumination func- 
tion describes the whole sphere being illuminated. This 
means that all allowed orbits lie in the disk and that 
the rest of the sphere does not contain broad line emit- 
ting gas. The second simulated dataset is also a thin 
disk viewed nearly face-on with a single value of energy 
and angular momentum, but for this case the illumina- 
tion function constrains the disk. We choose a value of 
#o close to 7r/2 so that orbits are allowed in the entire 
sphere. The third simulated dataset is a fully illuminated 
sphere with orbits that are also allowed in the entire 
sphere, so again 9q is close to tt/2. This is still an axisym- 
metric configuration, as the BLR gas density imposed by 
the kinematics depends upon the ^-coordinate. The true 
parameter values of the three simulated datasets used to 
test the kinematics model are shown in Table [3] 

We test each of the three simulated datasets assuming 
only one set of kinematics parameters. The parameter 
values inferred using our method are shown in Table 21 
while the full posterior pdfs are shown for all parameters 
of interest for the first simulated dataset in Figures [T2] 
and [T2J The posterior pdfs for the black hole mass, av- 
erage radius of BLR gas mass, and average width of the 
BLR gas mass are also shown for the second and third 
simulated datasets in Figures [T4l and [TBI They show that 
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Fig. 13. — Posterior pdfs for the first dynamical simulated 
dataset: face-on disk with the orbits confined to the disk. (Top) 
inclination angle, (middle) 8q, and (bottom) the joint pdf of 0q 
and the illumination angle. Notice in the joint pdf that 8o may 
only be larger than ~ 0.3 radians when the illumination angle is 
~ 0.3 radians, so the angular extent of the disk is well determined. 



0.16 



0.12 



CD 



0.08 



0.04 




6.95 7 7.05 7.1 7.15 
Log[Mbh/(solar masses)] 



SO. 15 



CD 




0.05 



14.03 14.06 14.09 
Log[Mean radius/m] 



0.16 



^ 0.12 



§ 0.08 



0.04 




13.55 13.6 13.65 
Log[Mean width/m] 



Fig. 14. — Posterior pdfs for the second dynamical simulated 
dataset: face-on disk with the orbits in the entire sphere. (Top) 
black hole mass, (middle) the average radius of the BLR gas mass, 
and (bottom) the average width of the BLR gas mass. 



the black hole mass, average radius, and average width 
of the BLR are well determined for all three simulated 
datasets. The angular parameters are also well deter- 
mined when physically possible. For example, for the 
first dynamics simulated dataset of a face-on disk with or- 
bits confined to the disk, the inclination angle and 0q are 
determined to within one or two grid points, while the il- 



lumination angle is only constrained to be > 0.3 radians. 
The illumination angle cannot be determined more ac- 
curately because the BLR gas emission only comes from 
the disk, so as long as the entire disk is illuminated the 
spectrum is not sensitive to further changes in the illu- 
mination angle. 

Finally, we also compute the velocity-resolved transfer 
functions for the three simulated datasets, shown in Fig- 
ure [THJ As expected, the transfer functions for the face- 
on disk configurations show little response at very small 
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Fig. 15. — Posterior pdfs for the third dynamical simulated 
datasct: sphere configuration with orbits allowed in the entire 
sphere. (Top) black hole mass, (middle) the average radius of the 
BLR gas mass, and (bottom) the average width of the BLR gas 
mass. 



lags, while the sphere configuration shows the highest 
intensity of response at small lags. The transfer func- 
tions for the face-on disk configurations are similar, but 
clearly lead to different line profiles, again illustrating 
the power of modeling the full dataset rather than just 
trying to model the transfer function. 

5. SUMMARY AND CONCLUSIONS 
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Fig. 16. — Velocity-resolved transfer functions for the three dy- 
namical simulated datasets: (top) face-on disk with orbits confined 
to the disk, (middle) face-on disk with orbits allowed in entire 
sphere, and (bottom) sphere configuration with orbits allowed in 
entire sphere. The red crosses show the response weighted mean 
lag in 10 velocity bins across the spectra. 



We introduce and test a new method for analyzing re- 
verberation mapping data of AGN by directly modeling 
the BLR. We illustrate our method by creating simple 
geometry and dynamical models of the BLR. Using a 
model of the BLR geometry to reproduce the integrated 
line flux timeseries from reverberation mapping data al- 
lows us to estimate the average radius of the BLR, as 
well as the mean width, illumination function, and incli- 
nation angle to the line of sight. Models of the BLR that 
include geometry and dynamical information allow us to 
additionally estimate the black hole mass and obtain an 
estimate of the extent to which the BLR gas orbits are 
confined to a disk or the whole sphere. 

Our method of analysis provides several advantages 
over previous methods. First, previous methods rely 
upon cross-correlation to obtain a mean radius for the 
BLR and a virial relation with unknown virial coefficient 
to obtain an estimate of the black hole mass. Our method 
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estimates the black hole mass self-consistently, without 
the need for a virial coefficient. Second, work model- 
ing reverberation mapping data has previously focused 
on modeling the velocity-resolved or unresolved transfer 
function. However the implications for the geometry and 
kinematics of the BLR are not clear for such analysis, as 
the transfer function is a function of the lag between 
the continuum and line emission. Instead of modeling 
the transfer function and then interpreting the transfer 
function in terms of a geometrical or dynamical model of 
the BLR, we focus on modeling the BLR directly. This 
allows us to extract more information and thus constrain 
the models more tightly. Finally, our fast method pro- 
vides estimates of the uncertainty in the model parameter 
values and can be used with numerical algorithms such 
as Nested Sampling that allow for model selection. Our 
main results can be summarized as follows: 

1. We create simulated datasets using the geometry 
model with known true parameter values and find 
that we can recover these values with uncertainties 
that depend upon the random uncertainty of the 
reverberation mapping data. We can recover the 
mean radius of the BLR to within ~ 0.1 dex and 
the mean width of the BLR to within ~ 0.2 dex for 
simulated data with an integrated line flux uncer- 
tainty of 1.5%. We can also place constraints on 
the inclination and illumination with uncertainties 
of ~ 0.2 radians for simulated data with face-on 
and spherical geometry configurations and 1.5% in- 
tegrated line flux uncertainty. Current integrated 
line flux uncertainties of about ~ 5% are on the 
edge of what would allow for successful recovery of 
more than just a mean radius for the BLR. 

2. We create simulated datasets using the dynamical 
model that consist of timeseries of a broad line pro- 
file and we compare them to mock spectra made 
using our model. Despite the larger number of free 
parameters in our dynamical model, we find that 
we can recover all the parameters physically possi- 
ble because the line profile is a stronger constraint 
on the model than the integrated line flux. We can 
recover the black hole mass and the mean radius of 
the BLR to within ~ 0.05 dex, for simulated data 
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with a line profile signal to noise ratio of ~ 4 per 
spectral pixel. We can also recover the mean width 
of the BLR to within ~ 0.1 dex and the inclination 
angle and illumination angle to within ~ 2 grid 
spacings over which the BLR density is defined. 

The small random uncertainties obtained in our tests 
of the simple geometry and dynamical models are partly 
due to the inherent assumption that our simulated data 
is drawn directly from the set of possible model con- 
figurations. In order to simulate the expected system- 
atic error in applying simple models to complicated real 
BLR systems, we have added substantial Gaussian noise 
to instances of the model in order to create our simu- 
lated datasets. The timeseries of line profiles, in the case 
of the dynamical model, is very constraining, and leads 
to the reduced random uncertainty in the mean radius 
and mean width of the BLR by a factor of two for the 
dynamical model, as compared to the geometry model. 
When applying the method to real data we expect larger 
uncertainties, owing to modelling errors. The uncertain- 
ties quoted here should therefore be considered as lower 
limits to the overall precision of the method for data of 
comparable quality. This emphasizes the importance of 
good quality data and increasingly more realistic models, 
for recovering detailed information about the BLR from 
reverberation mapping data. 

While we have created and tested both simple geome- 
try and dynamical models, our method is more general, 
allowing for use of any geometry or dynamical model that 
can be simply parameterized. We plan to expand our li- 
brary of models to include inflowing or outflowing BLR 
gas, which may be needed to explain some of the line pro- 
file asymmetries of current reverberation mapping data. 
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TABLE 1 

Simulated Geometry Data True Parameter Values 



Data 


Model 


Uncertainty* 


[10 14 m] 


<7 r 

[10 14 m] 


Inclination Angle 
[radians] 


Illumination Angle 
[radians] 


1 


Inclined Disk 


1.5% 


5 


0.3r mcan = 1.5 


0.79 


0.22 


2 


Inclined Disk 


5% 


5 


1.5 


0.79 


0.22 


3 


Edge-On Disk 


1.5% 


5 


1.5 


tt/2 


0.22 


1 


Face-On Disk 


1.5% 


5 


1.5 


0.0 


0.22 


5 


Shell 


1.5% 


5 


1.5 




tt/2 



Note. — Each simulated dataset consists of 60 line emission datapoints and the same 120 continuum emission 
datapoints, where the line emission timcscrics start half-way through the continuum timcscries. *Line flux uncertainty 
of the simulated dataset. 



TABLE 2 

Simulated Geometry Data Recovered Parameter Values 



Data 


Model 






Inclination Angle 


Illumination Angle 






[10 14 m] 


[10 14 m] 


[radians] 


[radians] 


1 


Inclined Disk 


4.53 ±0.47 


2.93 ±0.78 






2 


Inclined Disk 


4.81 ± 1.10 


2.24 ± 1.56 






3 


Edge-On Disk 


4.76 ± 0.65 


1.83 ±0.60 






1 


Face-On Disk 


4.93 ±0.67 


3.10 ± 1.31 


0.23 ± 0.13 


0.29 ±0.10 


5 


Shell 


4.48 ±0.52 


1.85 ±0.45 




1.22 ±0.29 



Note. — Results for 150,000 iterations using an MCMC algorithm. See Scction[4]for a discussion 
of why average values for the angular parameters are not quoted for most simulated geometry data- 
sets. 



TABLE 3 

Simulated Dynamics Data True Parameter Values 



Data 


Model 


<S/N>* 


M BH 
[1O 7 M ] 


Mean radius 
[10 14 m] 


Mean width 
[10 13 m] 


Inclination Angle 
[radians] 


Illumination Angle 
[radians] 


e a 

[radians] 


1 


Face-on Disk 


4.6 


1 


1.132 


4.484 


0.1 


tt/2 


0.3 


2 


Face-on Disk 


4.4 


1 


1.195 


4.022 


0.1 


0.3 


tt/2 


3 


Sphere 


3.5 


1 


1.195 


4.022 


0.1 


tt/2 


tt/2 



Note. — *Averagc signal to noise of the line flux profile. Each simulated dataset consists of 60 line emission profiles and the same 
120 continuum emission datapoints, where the line emission timeseries start half-way through the continuum timcscries. The simulated 
line emission profiles are created by taking the true model and adding gaussian noise with a variance of v — a. X Flux + j3, where 
a = 0.00018 and /3 = 0.025. 



TABLE 4 

Simulated Dynamics Data Recovered Parameter Values 



Data 


Model 


M BH 


Mean radius 


Mean width 


Inclination Angle 


Illumination Angle 


do 






[io 7 m ] 


[10 14 m] 


[10 13 m] 


[radians] 


[radians] 


[radians] 


1 


Face-on Disk 


0.95 ±0.05 


1.04 ±0.04 


4.27 ±0.05 


0.11 ±0.07 




0.31 ±0.02 


2 


Face-on Disk 


1.10 ±0.13 


1.12 ±0.04 


3.80 ±0.03 


0.12 ±0.06 


0.31 ± 0.02 


1.40 ±0.13 


3 


Sphere 


1.00 ±0.04 


1.21 ± 0.08 


3.95 ±0.04 


0.12 ±0.07 


1.46 ± 0.06 


1.50 ±0.05 



Note. — Results for 470 x 10 3 (data 1), 330 X 10 3 (data 2), and 110 x 10 3 (data 3) iterations using an MCMC algorithm. Sec 
Scction[4] for a discussion of why average values for the angular parameters are not quoted for the illumination angle of data 1. 



